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A new framework for evaluating hydrodynamic models of relativistic heavy ion collisions has been 
developed. This framework, a Comprehesive Heavy Ion Model Evaluation and Reporting Algo- 
rithm (CHIMERA) has been implemented by augmenting UVH 2+1D viscous hydrodynamic model 
with eccentricity fluctuations, pre-equilibrium flow, and the Ultra-relativistic Quantum Molecu- 
lar Dynamic (UrQMD) hadronic cascade. A range of initial temperatures and shear viscosity to 
entropy ratios were evaluated for four initial profiles, N part and N co u scaling with and without pre- 
equilibrium flow. The model results were compared to pion spectra, elliptic flow, and femtoscopic 
radii from 200 GeV Au+Au collisions for the 0-20% centrality range. Two sets of initial density 
profiles, Np ar t scaling with pre-equilibrium flow and N co u scaling without were shown to provide a 
consistent description of all three measurements. 



I. INTRODUCTION 

During the last decade a remarkable success has been 
achieved in the modeling and analysis of relativistic 
heavy ion collisions. Data from the Relativistic Heavy 
Ion Collider (RHIC) and more recently from the Large 
Hadron Collider (LHC) have surpassed the petabyte scale 
and include numerous classes of observables. The mod- 
els that have been most successful in describing the soft 
physics observables of particle spectra and collective flow 
measured at RHIC (see [THI] for a synopsis of the ini- 
tial results) are those that incorporate relativistic hydro- 
dynamics [5-7J coupled to a microscopic transport such 
as UrQMD [5J [5] ■ The relativistic hydrodynamic stage 
is used to describe the quark gluon plasma phase and 
its transition to hadronic matter, whereas the transport 
stage models the hadronic cascade which follows. Sub- 
sequent refinements to the hydrodynamic models to im- 
prove agreement with measured moments of the momen- 
tum anisotropy have incorporated shear viscosity [101 lllj 
and initial energy density fluctuations in the initial condi- 
tions [T2HTi] . The success of these models in reproducing 
the general features of particle spectra and elliptic flow 
has led to the conclusion that the QCD matter created 
in relativistic heavy ion collisions behaves very much like 
a fluid with small shear viscosity to entropy ratio. 

Although the comparisons between models and data 
are both interesting and compelling, there is not yet a 
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complete and rigorous quantification of the uncertainties 
in the parameters that would describe the quark gluon 
plasma created in heavy ion collisions at RHIC and the 
LHC. The greatest uncertainty lies in our understanding 
of the initial, pre-thermalized state of the collision. Dif- 
ferences in the initial density profile that alter the hydro- 
dynamic evolutions can lead to large uncertainties in the 
shear viscosity and initial temperature. The Equation of 
State (EoS) is calculated with increasing accuracy with 
lattice QCD (THIUS], but the interplay between the QCD 
EoS and the other parameters is yet be explored. Fur- 
thermore, most model comparisons are often limited to 
a subset of measurements that are deemed most relevant 
to a specific line of inquiry. Some success in uncertainty 
quantification has been achieved by studying the relation 
between two observables, such as the ratio of elliptic flow 
to eccentricity and the multiplicity density |17j . However, 
a complete determination of the properties of matter cre- 
ated in heavy ion collisions will require a comprehensive 
comparison between a hybrid hydrodynamic model and 
a set of physics observables. 

This paper describes a first step towards performing 
such a comparison, with the ultimate goal to fully con- 
strain the properties of the quark gluon plasma. We 
have developed a Comprehesive Heavy Ion Model Evalu- 
ation and Reporting Algorithm (CHIMERA) for system- 
atically comparing a set of hybrid hydrodynamic mod- 
els for heavy ion collisions spanning a range of initial 
parameters. This framework enables one to determine 
the optimal parameters and associated uncertainties that 
best describe a set of soft physics measurements, incor- 
porating both statistical and systematic errors. In the 
current implementation we use the 2D+1 viscous hydro 
code VH2 10J augmented with initial state eccentricity 
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fluctuations [T5] and pre-equilibrium flow [TS] to describe 
the hydrodynamic evolution, and the Ultra-relativistic 
Quantum Molecular Dynamics hadronic cascade code 
(UrQMD) [SI [S] to describe the hadronic transport. To 
compare to data we generate particle spectra and ellip- 
tic flow directly from the UrQMD output. We gener- 
ate femtoscopic correlation lengths, also referred to as 
HBT radii [201 [2T] , using the Correlation After-Burner 
(CRAB) code [23] • A chi-squared statistic is used 
to determine the best fit initial state parameters and as- 
sociated uncertainties for the measured results and er- 
rors. For this paper, we compare to published spec- 
tra, elliptic flow, and femtoscopic correlations for pions 
measured by the STAR and PHENIX collaborations for 
y/sNN=200 GeV Au+Au collisions in the 0-20% central- 
ity region. 

In Section [H] we provide a detailed description of the 
various components of the model, and in Section HI we 



explain the procedures used to generate the model results 
and to evaluate the agreement with the data. The results 
of these evaluations are reported in Section |IV| and po- 
tential systematic errors are discussed in Section [Vj In 
Section |VI| we give our conclusions and discuss future 
prospects and improvements. 



II. THE MODEL 
A. Initial Conditions and Hydrodynamic Evolution 

The primary component of the model we employ is the 
freely available 2+ ID viscous hydrodynamic code , VH2, 
developed by Luzum and Romatschke [TU]. VH2 numer- 
ically solves the Muller-Israel-Stewart equations with fi- 
nite shear viscosity in two dimensions assuming Bjorken 
boost invariant expansion along the longitudinal (beam) 
axis. We have modified the distributed version of VH2 in 
two ways: 1.) to account for the impact of initial density 
fluctuations on the average eccentricity for a given cen- 
trality, and 2.) to allow the initial state to be prepared 
with pre-equilibrium flow. 

The standard version of VH2 determines the initial 
density profile by calculating the Glauber overlap integral 
for the Wood-Saxon density distributions of the colliding 
nuclei. The initial conditions are therefore smooth, and 
the density profile may be set by participant or binary 
collisions scaling. We have modified the calculation of the 
initial conditions to account for initial eccentricity fluc- 
tuations, using the TGlauber Monte Carlo method [23] 
with the default parameters i?=6.38 fm, a=0.535, and a 
nucleon-nucleon cross-section cr/vjv=42 mb. Smooth dis- 
tributions were then obtained by averaging over 10,000 
events, in which each event was rotated to align the par- 
ticipant or binary-collision reaction plane along the X- 
axis. An event sample of this size will still produce fluc- 
tuations in the tails of the distribution that the VH2 code 
is not designed to handle. In order to remove these fluc- 
tuations the portion of the distribution below 3% of the 



maximum was fit to a 2-D Gaussian in X and Y, which 
was then used to replace the initial distribution. As be- 
fore, the initial density profile can be given by participant 
scaling, binary collisions scaling, or a superposition of the 
two. In this paper we use both participant and binary 
scaling in order to bracket the range of initial density dis- 
tributions commonly generated with a Glauber model. 

We have also introduced an option to prepare the ini- 
tial state with pre-equilibrium flow. Vredevoogd and 
Pratt have shown that for a system with a traceless stress 
energy tensor that obeys Bjorken boost-invariant scaling 
and for which the anisotropy in the spatial components is 
largely independent of the spatial coordinates, the flow 
can be expressed as a universal function of the energy 
component and time |19j . 
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As implemented in this work, the individual fluid cells 
(0.2 GeV -1 on a side) are given initial transverse ve- 
locities according to Eq. [I] after the initial energy den- 
sity distribution has been established. The addition of 
pre-equilibrium flow is expected to have the most impact 
on the elliptic flow and transverse radii. This can be 
seen in Fig. [TJ where the upper panels show the evolu- 
tion of the radial (left) and elliptic flow (right) through 
the mean transverse velocity and stress tensor asym- 
metry as a function of proper time. Radial and ellip- 
tic flow values are shown as a function of the proper 
time for start times of 0.6 fm/c (dotted) and 1.0 fm/c 
(solid). These are compared to VH2 run without pre- 
equilibrium flow for start times of 0.2 fm/c (double-dot- 
dashed), 0.6 fm/c (dot-dashed), and 1.0 fm/c (dashed). 
As the starting time is advanced, the systems that do 
not include pre-equilibrium flow begin to resemble the 
systems with pre-equilibrium flow, but they do not reach 
the same final values. Advancing the start time without 
pre-equilibrium flow also advances the freeze-out time by 
a similar amount. Note that the two systems prepared 
with pre-equilibrium flow do not show much dependence 
on start time, a result that is consistent with the premise 
of a universal pre-equilibrium flow. The lower panels of 
Fig. [I] show the mean radial position and radial asymme- 
try for the same conditions. Again the earlier start times 
without pre-equilibrium flow approach, but do not match 
the values associated with pre-equilibrium flow, which 
show little dependence on start time. For the compar- 
isons to experimental data that follow, we adopt a start 
time of 1.0 fm/c with and without inclusion of the pre- 
equilibrium flow. This start time is consistent with the 
choice in [TU] for Glauber initial conditions. 

The hydrodynamic evolution then proceeds as with 
the standard version of the VH2 code, using the QCD- 
inspired Equation of State based upon the work of 
Laine and Schroeder [25] that interpolates between the 
hadronic resonance model and the perturbative calcula- 
tion. It provides a reasonable albeit imperfect approx- 
imation to the cross-over transition that is now calcu- 
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FIG. 1. The proper time evolution of mean radial velocity (upper left) and stress tensor asymmetry (upper right) flow for 
VH2 2+1D hydrodynamic evolution for Au+Au collisions with participant Glauber initial conditions for a b = 4.4 fm impact 
parameter and corresponding b = initial temperature of 300 MeV. Contours are drawn for initial proper start times, r=0.2, 
0.6, and 1.0 fm/c without pre-equilibrium flow and for r=0.6 and 1.0 with pre-equilibrium flow as indicated in the legend by 
the ft symbol. The mean transverse radial position (lower left) and eccentricity (lower right) are also shown. Note that all 
means are energy weighted. 



lated non-perturbatively using lattice QCD. A thorough 
study of the implications of various lattice calculations 
and their equation of state parameterizations is left for 
future investigation. 



B. Freeze-out and Hadronic Cascade 

The hydrodynamic evolution in VH2 ceases when the 
temperature of a given cell falls below a specified freeze- 
out value. The final energy densities are converted to 
final state particles following the prescription of Cooper 
and Frye |26| with corrections for the shear viscosity im- 
plemented according to the method developed by Pratt 
and Torrieri |27) . In this method, the particle number 
and momentum distributions are determined by Monte 
Carlo sampling for non-viscous equilibrium distributions, 
and the momenta are subsequently rescaled according to 
Eq.0 



Here A^- is proportional to tt^ , the Israel-Stewart correc- 
tion to the stress energy tensor. The constant of propor- 
tionality is chosen to reproduce the second order viscous 
corrections to the final particle distributions for small 
values of ir. The list of final state particles is selected to 
match the complete set of known particles used by the 
UrQMD code. For each CHIMERA parameter set a file 
of 5000 events is generated in the OSCAR-97 format [28] . 
For comparisons to elliptic flow and radii, an additional 
15000 particles are generated in order to achieve smaller 
statistical errors in the model results prior to fitting. We 
use a switching temperature of T sw =165 MeV to end the 
hydrodynamic evolution and generate the particles for 
the subsequent hadronic cascade stage of the model. As 
in [IT], the value of T sw was chosen to be as close as 
possible to, but still larger than the transition tempera- 
ture. The freeze-out particle times are converted to for- 
mation times and the particles are back-propagated to 
zero time using the standard oscar2u program provided 
by the UrQMD developers upon request. 



V 3 -> Pi + hjPj- 



(2) 



For this paper we use UrQMD v2.3 to model par- 
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ticle interactions that follow the hydrodynamic freeze- 
out stage. UrQMD is run with default switches except 
for modifications needed to read the OSCAR-97 format- 
ted input file. We also disable the unstable particle 
decay after the final time step, in order to match the 
full set of particle spectra measured by the experiments. 
The CHIMERA framework can also be used without the 
hadronic cascade. In this case, the list of final state par- 
ticles is equal to the set of stable particles from the 2008 
listing of the Particle Data Group [29] and a lower freeze- 
out temperature T stu =140 MeV is used. 



C. Post Processing 

The final state particle distributions from UrQMD are 
then used to construct the observables that can be di- 
rectly compared to published experimental data. In this 
work we restrict our comparison to invariant particle 
spectra, the second coefficient of the transverse momen- 
tum anisotropy with respect to the event reaction plane 
i>2, and the Bertsch-Pratt femtoscopic radii, i?i on g, -Rsidc, 
R ou t- The transverse momentum spectra are calculated 
as invariant cross-sections using bins of 0.1 GeV/c in 
the range 0.2-1.5 GeV/c within the rapidity interval of 
\y\ < 0.5. The elliptic flow is calculated within the same 
transverse momentum region. We use CRAB to generate 
the three-dimensional femtoscopic correlation functions, 
including the strong interaction and statistical interfer- 
ence. The Coulomb interaction is neglected so that a 3D 
Gaussian can be used to fit the correlation function to 
generate the radii. The correlation functions are binned 
in 0.2 GeV intervals. The fits are performed directly 
on the correlation weights, to avoid the time and cpu- 
consuming process of constructing an event mixed back- 
ground for each momentum bin. 

The cumulative effect of each modification to the VH2 
code is shown in Figures [2] and [3j Fig. [2] shows the 
charged particle yields and mean transverse momentum 
((pr)) for pions (filled triangles), kaons (filled inverse 
triangles), and protons (open triangles). The unmod- 
ified VH2 results are shown as solid circles with con- 
necting lines, and correspond to a binary scaling optical 
Glauber model input for an initial central temperature 
of T cent — 333 MeV, r]/s = 0.08 and a freeze-out tem- 
perature of 140 MeV. These parameters were chosen to 
match the upper panels of Fig. 7 of [3D]. Here we fol- 
low the VH2 convention in reporting the equivalent tem- 
perature, T centl for a central collision (zero impact pa- 
rameter). As in |30j . the model results are compared to 
PHENIX measurements from [3T] . The results in Fig. [2] 
make use a different freeze-out routine, described above, 
that was written to accommodate the need to write OS- 
CAR formatted events for input to CRAB and UrQMD. 
Therefore the results may differ slightly from the pub- 
lished VH2 results. This figure shows the cumulative 
effect of introducing eccentricity fluctuations (open cir- 
cles), pre-equilibrium flow (open squares) and UrQMD 



with a switching temperature of T sw — 165 MeV (filled 
squares). The incorporation of the eccentricity fluctua- 
tions has only small impact on the particle yields and 
(pr), but the addition of pre-equilibrium raises the (pr) 
by as much as 20%. The hadronic cascade phase leads 
to increased proton yields and decreased (pt) for pions 
and kaons; it is likely that both result from additional 
collisions in UrQMD involving baryonic resonances. In 
adding features to the model, we make no attempt to 
modify the initial parameters to achieve better agree- 
ment with the data. This will be the main focus of the 
results section. 

Similar comparisons for the elliptic flow for charged 
hadrons and sidewards and outwards radii for pions are 
shown in Fig. [3] Each modification to VH2 serves to in- 
crease the V2- For N part > 200 the model results come 
closer to the comparison data measured by STAR [3J, 
however for more peripheral collisions the addition of pre- 
cqiuilibrium flow causes the value of i>2 to overshoot the 
data by a significant margin. For very peripheral colli- 
sions the anisotropy of the stress energy tensor may no 
longer be independent of the spatial coordinates, thereby 
violating the conditions required for universal flow [19j . 
The hadronic cascade phase leads to a reduction in the 
value of vi , when compared to a pure hydrodynamic evo- 
lution with lower freeze-out temperature. 

The right panel of Fig. [3] shows the cumulative effects 
on the transverse radii, R ou t and R S ide- Note that each 
modification leads to a reduction in R ou t and an increase 
in R s idei bringing the two closer to each other and to 
the ratio observed in the data. The measurements by 
STAR and PHENIX for 200 GcV Au+Au are omitted 
from this figure for clarity in order to better reveal the 
trends, which are consistent with the one-dimensional rel- 
ativistic viscous hydro results documented in |32j . 



III. MODEL EVALULATION 

A. Model Parameterizations 

Because the published experimental results do not gen- 
erally conform to a uniform binning, the model results 
are fit to a functional form prior to evaluating the overall 
agreement with the data. This step also simplifies the 
treatment of systematic errors in the data. We selected 
a purely empirical fitting function in order to avoid an 
bias in treating the model results prior to comparison to 
experiment. The pt dependence of i>2 and the fcy depen- 
dence of the femtoscopic radii are both well described by 
a Chebyshev polynomial. However, polynomials in gen- 
eral are not easily adapted to the large dynamic range 
characteristic of spectra. The spectra require an expo- 
nential form to describe the data, and we achieved satis- 
factory agreement by multiplying a fifth order Chebyshev 
polynomial by an exponential in transverse mass that has 
been shown to fit a wide range of particle spectra |31) . 
The functional form use to fit the model spectra is given 
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FIG. 2. dN/dy (left) and mean transverse momenta (right) for pions (red), kaons (green), and protons (blue) for the unmodified 
VH2 (solid connected circles), and for the cumulative effect of adding eccentricity fluctuations (open circles), pre-equilibrium flow 
(open squares) and switching to the UrQMD hadronic cascade (solid squares) at T sw — 165 MeV. Model results are compared 
to PHENIX measurements for pions (filled triangles), kaons (inverted filled triangles), and protons (open triangles). In this 
figure the open-circle symbols for the eccentricity fluctuations are often obscured by filled-circle results from the unmodified 
VH2 results. 




FIG. 3. Mean elliptic flow for VH2, adding eccentricity fluctuations, pre-equilibrium flow, and the UrQMD cascade after-burner 
(left). R a ide and R ou t femtoscopic radii for the same conditions (right). 



byEq.^ 



B. Data Comparisons 



v(Pt) = ^2aiTi(p T ) 



exp [-(iB T - m )/T] 
2irT(T + m ) 



(3) 



We count on the Chebyshev multilpier to cancel any bias 
that may result from incorporating the slope parameter, 
T, into the fit. We have verified that these functional 
forms provide a good description of the model results, 
with Xdof c l° se to unity for most systems. 



To evaluate the model parameters we calculate a chi- 
squared statistic for each combination of model parame- 
terization and data set. We have selected four different 
initial states for evaluation, including both participant 
(Npart) and binary collision (N co u) scaled energy density 
profiles, each prepared with and without pre-equilibrium 
flow. The initial Glauber distributions were generated for 
a fixed impact parameter of 4.4 fm, corresponding to the 
0-20% centrality bin for 200 GeV Au+Au collisions. The 
initial viscosity to entropy ratios were sampled in steps 
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of 0.08, from a lower bound of 0.001 (a non-zero value 
is required for numerical stability) to an upper bound of 
0.48. The temperature was varied in steps of 5 MeV for 
spectra and 20 MeV for flow and radii, with the ranges 
adjusted iteratively to reach above and below the optimal 
parameter values for each set of initial conditions. The 
full set of input parameters and data comparison sets are 
listed in Table [I] A moderate centrality range of 0-20% 



this assumption, the modified \ 2 formula reduces to, 



Initial Profile 

Np ar t 


T)/s 

0.001-0.48 


ry-ispectra 
cent 

280-325 MeV 


rpv 2 — f emt o 
cent 

260-380 MeV 


Npart pre-equ. flow 


0.001-0.48 


270-315 MeV 


260-380 MeV 


N con 


0.001-0.48 


330-370 MeV 


300-420 MeV 


N coH pre-equ. flow 


0.001-0.48 


310-350 MeV 


300-420 MeV 



TABLE I. Initial VH2 density profile and range of central 
temperatures (T cen t) and T]/s values. All runs used a fixed 
4.4 fm impact parameter with (iVp ar t}=276, t—1.0 fm/c, the 
default QCD-inspired EOS, and T s „=165 MeV. 



was selected to avoid large corrections to the v 2 associ- 
ated with non-flow effects in the data [33] and to match 
existing pion femtoscopy, v 2 , and spectra. The pion spec- 
tra are compared to the 0-20% centrality measurements 
by PHENIX [3JJ and STAR [3J. CHIMERA evaluates 
models to data by centrality using (N part ). For this 
impact parameter we calculate a value of (N part )=276, 
whereas PHENIX reports a value of 286 and STAR re- 
ports a value of 282 for their respective 0-20% centrality 
ranges. The nearest pion radii for PHENIX are for the 
centrality range of 0-30%, but the presence of two pi- 
ons in one of the PHENIX central arms imposes a slight 
centrality bias and (N part }=281 for these data [3S]. The 
nearest centrality bin for the STAR pion radii is 5-10%, 
which have (N part )—298 [36]. Based on the N part scaling 
analysis in [37], this leads to a difference in radii of less 
than 2%, well below the errors of cither measurement. 
For the v 2 we compare to two sets of measurements by 
PHENIX in this centrality range, a combined analysis of 
pions and kaons |38] , and a recent analysis of pion v 2 at 
higher px [3j?j. The STAR v 2 results for pions are also 
from the 5-10% centrality range [3J>] . From the centrality 
dependence shown in Fig. [3j a reduction of 26 in (N part ) 
can produce an increase in v 2 ~8%, which is compara- 
ble to the reported errors and may be large enough to 
influence the results. The impact of this discrepancy is 
examined further in Section IVl 

To incorporate the systematic errors when evaluating 
the chi-squared statistic, we follow the procedure defined 
in 40J. Because it is generally not feasible to calculate 
a full co-variance matrix for the systematic errors, the 
authors of [ID] make a set of simplifying assumptions that 
correspond to different types of systematic errors. For 
this analysis we assume all systematic errors assumed to 
be of type B, in which the systematic errors are assumed 
to be fully correlated in the fractional deviation. With 



E 



of (1 + e b a b Jyi) 2 



(4) 



where a bi represents the systematic error, and e b repre- 
sents the fraction of the error by which the set of cor- 
related measurements move in tandem. The full proce- 
dure for evaluating the modified x 2 and the derivation 
are given in the appendix of [40 . 



IV. RESULTS 

The full set of CHIMERA model evaluation results are 
shown in Figures [4 ■ 15 



For each of these figures, the 
leftmost panel shows the evaluation with fixed r\js set 
to the conjectured minimum value of 0.08. The middle 
panel shows the variation with eta/ s for fixed tempera- 
ture, set to the approximate minimum. The Chebyshev 
fits to the model results are shown and the experimen- 
tal data are plotted with both statistical error bars and 
systematic error boxes. The calculated values of x\ f 
are reported separately for each experiment. Because 
the wide dynamic range of the spectra can obscure sub- 
tle differences, the ratios of data to model fits are also 
shown. The upper rightmost panel shows the sum of the 
X 2 over all data sets divided by the total degrees of free- 
dom (dof). A paraboloid fit to the x\of distribution is 
shown immediately below each distribution. These fits to 
the Xdof distributions will be used to extract confidence 
contours for each measurement. 

Fig. HI shows the model fits, data, and Xdof distri- 
bution Tor the pion spectra from N part scaling without 
pre-equilibrium flow. In general, the calculated spectra 
are harder (less steep) than the data, a condition that 
will persist for the other initial conditions. The spec- 
tra harden further with increasing rj/s, a feature evident 
in [10] but not often shown due to the propensity to use 
the spectra to constrain only the temperature. However, 
this phenomena leads to an inverse correlation between 
Xdof mm i ma i n T and 77/3, as shown in the right panel 

of Fig. g 

In Fig. [5] one can see that the addition of pre- 
equilibrium flow serves to further harden the spectra. 
Some flexibility to find a reasonable x\ f minimum is 
derived from the freedom for the lower pt spectra from 
STAR to move independently from the PHENIX data, 
each within their systematic errors. However, the curva- 
ture for the Xdof distribution is noticeably steeper. The 
parameters for the paraboloid fits to these and all other 
measurements are listed in Table UTI 

The model comparisons for the N co u scaling are shown 
in Fig. [6] The N co u scaling results require a higher 
temperature to match the pion multiplicity, which fur- 
ther exacerbates the difference in slope. When the pre- 
equilibrium flow is added in Fig. [7] the agreement with 
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FIG. 4. Model evaluation of pion spectra with N par t scaling for fixed T (left), fixed r)js (middle), with total Xdof distribution 
(upper right) with paraboloid fit (lower right). 




FIG. 5. Model evaluation of pion spectra with N par t scaling and pre-equilibrium flow for fixed T (left), fixed r]/s (middle), 
with total Xdof distribution (upper right) with paraboloid fit (lower right). 



data is poorer. The Xdof values f° r this minima are sig- 
nificantly higher and the distribution is quite steep. 

Figures [8] and [9] show the model evaluations for the 
elliptic flow, with the absence and addition of pre- 
equilibrium flow, respectively. Comparisons are shown 
separately for pions from STAR and PHENIX, and ear- 
lier results from PHENIX that combine pions and kaons. 
Without the addition of pre-cquilibrium flow, the min- 
imum occurs near zero rj/s and at higher values of the 
temperature. With the addition of pre-equilibrium flow 
higher values of 77/s and lower temperatures can be ac- 
commodated. In contrast to the model comparisons for 
spectra, the x\of distributions for v-i evaluations exhibit 
more fluctuations, an indication that the event sample 
size should be increased beyond 20,000 events for future 
evaluations. 

The elliptic flow comparisons for N co u scaling are 
shown in Fig. |10| The larger pressure gradients of the 
N co u scaling push the minima to higher rj/s values, into 
a region of I < Anrj/s < 2, consistent with the initial 
VH2 results using Glauber N co u geometry [10] and the 
recent work of Song et al. [T7]. The addition of pre- 



set of concentric ellipses labeled corresponding min imum 

shows 



equilibrium flow shown in Fig. II pushes this minimum 
even higher. However, the Xdof distributions are quite 
broad. Aside from the exclusion of the high T, rj/s re- 
gion without pre-equilibrium flow and the low rj/s region 
when pre-equilibrium flow is added, there appears to be 
limited constraining power from comparing to the pt de- 
pendence of one particle species at a single centrality. 

In Fig. [12] we show for the first time a direct compar- 
ison of VH2+UrQMD values of Ri ong , Rside, an d R ou t 
compared to experimental data for 200 GeV Au+Au col- 
lisions. The comparison displays the characteristic dis- 
crepancy (referred to as the HBT puzzle [41] ) in which 
Rout trends above the data and R S ide trends below. Bet- 
ter agreement is obtained for lower temperatures, where 
a minimum value of Xdof below 20 is achieved in the 
T=260 MeV bin. However, this is also below the tem- 
perature constraints set by the spectra. When pre- 



equilibrium flow is added, as shown in Fig. 13 the overall 



agreement improves significantly, and the Xdof minimum 
occurs at temperature of 320 MeV, within the range of 
values allowed by the spectra. The femtoscopic radii have 
little to no constraining power on the shear viscosity in 
this evaluation. 

The model evaluation for N co u scaling pion radii re- 
sults shown in Fig. 14 are similar in feature to the N part 



scaling results. The Xdof minimum is shifted to higher 
temperature, but by an amount similar to the magni- 
tude of the shift that occurs for the spectra. Fig. 15 



shows the agreement improving again with the addition 
of pre-equilibrium flow. 

To examine the relationship between the Xdof distri- 
butions for each measurement more closely we show in 



Fig. 16 the one- and two-sigma contours obtained from 
the paraboloid fits. The two-sigma contours are solid and 
one-sigma contours are dot-dashed. Contours are drawn 
for the femtoscopic radii, elliptic flow, spectra, with each 
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value of Xdof achieved for each evaluation. Fig. 
the evaluations for N part scaling (left) and iT^li scal- 
ing (right) with pre-equilibrium flow (top) and without 
(bottom) . The regions that fall within the two-sigma con- 
tours for all measurements are shaded. This occurs for 
N par t scaling with pre-equilibrium flow as well as for the 
N co u scaling without. These are also initial conditions 
that lead the best overall agreement with the three data 
sets. The N part evaluation in the lower left quadrant of 
Fig. |16| provides a different perspective on what has been 
referred to as the "HBT puzzle" [? ]. The optimal tem- 
peratures for radii, spectra, and flow are inconsistent. 
Selecting the "best fit" temperature by eye leads to a 
poorer agreement with the radii as well as the flow. The 
addition of pre-equilibrium flow shifts the HBT minimum 
to higher temperature, and broadens the v2 distribution 
to achieve a consistent, albeit still imperfect agreement 
between the different measurements. The use of N co u 
scaling has a similar effect, but adding pre-equilibrium 
flow to this scaling overcompensates. 

Although the existence of a two-sigma overlap region is 
encouraging, it is obvious that we have not yet achieved 
a truly acceptable evaluation of x 2 dof ~ 1 for a single 
measurement, and even if this were achieved, it would 
be necessary to evaluate a greater breadth of measure- 
ments and centrality bins make a compelling case for 
the validity of the model and initial conditions. How- 
ever, when this is achieved, the maximum constraining 
power will be realized by combining the chi 2 evaluations 
from all measurements. Fig. [17] shows the distributions 
for the sum over x 2 f° r each of the three evaluations: 
spectra, flow, and radii, divided by the total degrees of 
freedom. A close inspection reveals that these distribu- 
tions are nearly identical to the Xdof distributions for the 



spectra show in Figures |4]|7j As expected, the x sums 
are completely dominated by the spectra. 

To achieve a better balance between the three mea- 
surements, we repeat the x 2 sum with the spectra given 
a 10% weight relative to the other measurements. These 
distributions are shown in Fig. 18 These distributions 



are still quite close to the spectra, the parameters listed in 
Table [TT] are slightly different, indicating that information 
from the elliptic flow and femtoscopic radii evaluation is 
being given greater weight. The overall chi 2 do f for the 
weighted sum, denoted as Suniio in the table, is below 
20 for the conditions without pre-equilibrium flow, and 
less than 10 for the initial conditions that include it. A 
comparison of the one- and two-sigma contours for the 
weighted sum x\of distributions is shown in Fig. 
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V. SYSTEMATIC ERRORS 

The incorporation of systematic errors from the mea- 
surements is an important component of the evaluation 
that has been presented. A separate source of systematic 
errors comes from the models employed. These include 
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FIG. 7. Model evaluation of pion spectra with N co u scaling with pre-equilibrium flow for fixed T (left), fixed r]/s (middle), 
with total Xdof distribution (upper right) with paraboloid fit (lower right). 
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FIG. 8. Model evaluation of pion elliptic flow with N par t scaling for fixed T (left), fixed r//s (middle), with total Xdof distribution 
(upper right) with paraboloid fit (lower right). 




FIG. 9. Model evaluation of pion elliptic flow with N par t scaling with pre-equilibrium flow for fixed T (left), fixed rj/s (middle), 
with total Xdof distribution (upper right) with paraboloid fit (lower right). 
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FIG. 10. Model evaluation of pion elliptic flow with Ncoll scaling for fixed T (left), fixed r//s (middle), with total Xdof 
distribution (upper right) with paraboloid fit (lower right). 




FIG. 11. Model evaluation of pion elliptic flow with N co u scaling with pre-equilibrium flow for fixed T (left), fixed rj/s (middle), 
with total Xdof distribution (upper right) with paraboloid fit (lower right). 
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FIG. 12. Model evaluation of pion radii with N par t scaling for fixed T (left), fixed r]/s (middle), with total Xdof distribution 
(upper right) with paraboloid fit (lower right). 




FIG. 13. Model evaluation of pion radii with N par t scaling with pre-equilibrium flow for fixed T (left), fixed rj/s (middle), with 
total Xdof distribution (upper right) with paraboloid fit (lower right). 
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FIG. 15. Model evaluation of pion radii with N co u scaling with pre-equilibrium flow for fixed T (left), fixed r)/s (middle), with 
total Xdof distribution (upper right) with paraboloid fit (lower right). 
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Initial profile 


Evaluation 


Major axis 


Minor axis 


angle (deg) 


T ■ 

J- mm 


V/ S min 


X 2 /d0f m in 


N part 


HBT 


0.211 


0.0114 


0.28 


0.266 


0.5 


13.9 


v 2 


0.088 


0.0134 


6.52 


0.347 


0.0 


11.5 


Spectra 


0.153 


0.0080 


2.3 


0.297 


0.29 


2.7 


Sum 


0.126 


0.0049 


2.18 


0.297 


0.22 


12.7 


Sum 10 


0.127 


0.0074 


1.24 


0.29 


0.21 


19.2 


N con 


HBT 


0.237 


0.0121 


1.01 


0.317 


0.52 


14.1 


v 2 


279 


0.0318 


19.2 


0.324 


0.12 


3.9 


Spectra 


210 


0.0023 


2.27 


0.348 


0.0 


4.4 


Sum 


620 


0.0032 


2.15 


0.341 


0.18 


10.6 


Sumio 


13.8 


0.0068 


2.01 


0.334 


0.28 


14.4 


prpre- flow 
part 


HBT 


0.374 


0.0105 


2.51 


0.273 


0.52 


4.8 


v 2 


0.335 


0.039 


14.2 


0.322 


0.0 


4.8 


Spectra 


0.094 


0.0019 


3.11 


0.303 


0.0 


4.9 


Sum 


0.128 


0.0027 


2.96 


0.302 


0.0 


6.3 


Sumio 


0.246 


0.0059 


2.69 


0.300 


0.0 


6.4 


A/pre- flow 
iV co(( 


HBT 


0.434 


0.013 


2.52 


0.361 


0.24 


4.3 


v 2 


0.156 


0.0552 


19.3 


0.29 


0.30 


3.5 


Spectra 


0.086 


0.0015 


3.05 


0.341 


0.0 


29.8 


Sum 


0.138 


0.0022 


3.05 


0.341 


0.0 


18.0 


Sumio 


0.182 


0.0050 


2.91 


0.328 


0.26 


8.4 



TABLE II. Ellipse parameters determined from paraboloid fits to Xdof distributions in T and r]/s for each set of initial condition 
profile and for each measurement. 



the systematic errors of Glauber model used to generate 
the initial conditions, the parameterization of the initial 
state flow, the choice of hydrodynamic solvers, the freeze- 
out conditions, and cascade model assumptions. A full 
investigation of these effects is well beyond the scope of 
the current evaluation and is left for future study. How- 
ever, there are additional set of systematic errors that 
are specific to the CHIMERA framework: 

• the criteria for centrality matching 

• the use of Chebyshev polynomials to fit momentum 
dependence of model results, 

• the momentum range over which evaluations were 
performed, 

• the paraboloid functions fit to the Xdof distribu- 
tions in T and 77/ s, 

• the range in T and r\j s used for the paraboloid fits. 

These items are also left for future study, except for the 
centrality matching criteria. As mentioned earlier, there 
difference between the model centrality of < N part =276 
and < Np ar t >=298 for the femtoscopic radii and elliptic 
flow measurements by STAR is large enough to have a 
significant impact on the evaluation procedure. To gauge 
the impact of this effect on the Xdof evaluations of fem- 
toscopic radii we repeat the analysis after performing a 
linear interpolation of the radii and V2 measurements to 



< N part =276. For the radii, the interpolation is a func- 
tion of Npar t using parameters derived from fits to the 
the centrality dependence |37j . For the v 2 measurements, 
the interpolation was performed using the integrated V2 
centrality dependence measured in [36] . As stated previ- 
ously, these adjustment amount to a 2% decrease for the 
STAR radii and an 8% increase in v 2 . The CHIMERA 
evaluation is repeated for only the N part scaling with 
pre-equilibirum flow, which was shown to have the best 
overall agreement with the data. The net effect of apply- 
ing this centrality adjustment is shown in Fig. [20] and the 



parameters for the paraboloid fits are listed in Table III 



Evaluation 


Major 


Minor 


angle 


Txain 


??/Smin 


Xdof 


HBT°* 


0.366 


0.0103 


2.52 


0.265 


0.52 


4.8 


HBT 


0.374 


0.0105 


2.51 


0.273 


0.52 


4.8 


yadj 


0.180 


0.0315 


15.5 


0.322 


0.0 


6.3 


v 2 


0.335 


0.0393 


14.2 


0.322 


0.0 


4.8 


adj 

Sum 10 J 


0.197 


0.0059 


2.69 


0.298 


0.0 


7.4 


Sumio 


0.246 


0.0059 


2.69 


0.300 


0.0 


6.4 



TABLE III. Ellipse parameters determined from paraboloid 
fits to Xdof distributions in T and r]/s for N par t scaling with 
pre-equilibrium flow and with the centrality adjustment de- 
scribed in the text. 

The differences between the unadjusted and adjusted 
centrality evaluations are minor. The best fit tempera- 




FIG. 16. Xdof contours for pions spectra, elliptic flow, radii, and a weighted sum in which the spectra Xdof are weighted by 
10% relative to the other measurements. The four panels show contours for N par t scaling (upper left), N co ii scaling (lower left), 
Np ar t with pre-equilibrium flow (upper right), and N co u with pre-equilibrium flow (lower right). The minimum Xdof values are 
labeled for each fit, and regions which overlap at the two-sigma level are shaded. 



ture for the radii shifts down by 8 MeV, and the contour 
for the V2 becomes narrower by 40% along the major axis 
and 20% along the minor axis. The overall Xdof mm i- 
mum rises slightly. We conclude that the overall central- 
ity match between the model and data is sufficient for 
this limited evaluation, however, this issue will require 
further scrutiny as the CHIMERA evaluations become 
more precise. 



VI. CONCLUSIONS 



We have performed a simultaneous \ 2 evaluation of 
three measurements: spectra, elliptic flow, and femto- 
scopic radii from 200 GeV Au+Au collisions for a the 0- 
20% centrality bin, using using an augmented version of 
the VH2 viscous 2D+1 hydrodynamic model that incor- 
porates pre-equilibrium flow and initial state eccentricity 
coupled to the UrQMD hadronic cascade. The evaluation 
was performed for four sets of initial conditions: N part 
and N co u initial density profiles with and without the 
addition of pre-equilibrium flow. The x 2 evaluation was 
performed for measurements by the PHENIX and STAR 
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FIG. 17. Equal weight sum of Xdof f° r pion spectra, elliptic flow, and femtoscopic radii. From left to right the distributions 
are for N par t scaling, N co u scaling, N par t scaling with pre-equilibrium flow, and N co u scaling and pre-equilibrium flow. The 
summed Xdof distributions are on top, and the paraboloid fits are below. 




FIG. 18. Sum of Xdof f° r N part scaling with spectra given 10% weight. From left to right the distributions are for N par t 
scaling, N co u scaling, N part scaling with pre-equilibrium flow, and N co u scaling and pre-equilibrium flow. The weighted sum 
Xdof distributions are on top, and the paraboloid fits are below. 



collaborations, and include both statistical and system- 
atic errors. For two of the initial conditions, N part with- 
out pre-equilibrium flow, and N co u with pre-equilibrium 
flow, the constrained regions of T and rj/s were mutu- 
ally exclusive. However, for the N part scaling with pre- 
equilibrium flow and N co u without, the constraints re- 



gions overlap at the level of 2<r. This is also the first 
time that pion radius measurements have been success- 
fully reproduced by a 2D+1 viscous hydrodynamic model 
coupled to a hadronic cascade. 

We regard the main conclusion of this work to be 
the successful demonstration of an evaluation technique 
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T (GeV) 

FIG. 19. Xdof contours for the weighted sum in which the 
spectra are given a 10% weight relative the x 2 distributions 
from elliptic flow and femtoscopic radii. One- and two-sigma 
contours are drawn with solid lines for the N par t (left) and 
Ncoii scaling, and with grey dashed lines for the addition of 
pre-equilibrium flow. 




T (GeV) 

FIG. 20. Xdof contours for pions spectra, elliptic flow, radii, 
and a weighted sum for N par t scaling with pre-equilbrium 
flow with an additional adjustment to improve the centrality 
match with the STAR femtoscopic radii and elliptic flow. 



The centrality matching has room for improvement, the 
equation of state is not yet specified according to the 
most recent lattice QCD results, and the full range of 
initial state profiles has not been explored. There are 
also additional model parameters, such as T star t and T sw 
that have not been explored. In addition, the compar- 
isons to experimental data were limited, covering only a 
single particle species, and a single centrality range. We 
expect to address these issues in the near future. 

Perhaps a more severe limitation on the current imple- 
mentation of CHIMERA is the inability of VH2 to handle 
non-smooth initial conditions. Qiu and Heinz have shown 
that [H] a proper accounting of event by event fluctua- 
tions is required to compare to non-central (>60%) col- 
lisions is also needed to generate the higher order har- 
monics that promise to improve our ability to constrain 
r)/s [42j[43]. In order for CHIMERA or a similar frame- 
work to ultimately succeed, a more sophisticated 3D+1 
hydrodynamic code will need to be used, and its perfor- 
mance will need to be evaluated against a wide range of 
physics observables. 
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for constraining the multi-parameter space of heavy ion 
models by comparing to multiple data sets. Working with 
a relatively small data sample, we were able to achieve 
a set of constraints for two sets of initial conditions, and 
to exclude two others. However, the current implemen- 
tation of CHIMERA should be regarded as incomplete. 
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